Update

We will update this tutorial when necessary. Readers can access the latest version in our GitHub repository.

If you have any questions, errors or bug reports, please contact Pietro Pollo (pietro_pollo@hotmail.com) or Shinichi Nakagawa (snakagaw@ualberta.ca).

Introduction

This online material is a supplement to our paper “Beyond sex differences in mean: meta-analysis of differences in skewness, kurtosis, and correlation”. You will see how to calculate the new effect size statistics we have proposed and how to use them in a meta-analytical model using the metafor package in R.

Content

In this online material, we provide details on the simulations we ran to evaluate the effectiveness of our proposed effect size statistics (and their associated sampling error). We also show how to (1) calculate our newly proposed effect sizes (\(\Delta sk\), \(\Delta ku\), \(\Delta Zr\)) and (2) exemplify their use with data from the International Mouse Phenotyping Consortium.

Simulations

To understand the effectiveness of our proposed effect size statistics at capturing differences in skewness, kurtosis, and correlation between two groups we ran a series of simulations. We simulated two distinct groups that varied in distributional moments (i.e., \(\Delta sk\) and \(\Delta ku\)) or in the correlation within each group (i.e., \(\Delta Zr\)). We used the Pearson function from the moments (Komsta and Novomestky 2022) package to control aspects of skewness and kurtosis in each group. For simulating differences in correlations between groups we used the mvnorm function in the MASS package [].

Prerequisites

Loading packages

Our tutorial uses R statistical software and existing R packages, which you will first need to download and install.

If the packages are archived in CRAN, use install.packages() to install them. For example, to install the metafor , you can execute install.packages("metafor") in the console (bottom left pane of R Studio).

Version information of each package is listed at the end of this tutorial.

if (!require("pacman")) {install.packages("pacman")}
pacman::p_load(corrr,
               DT,
               ggdist,
               ggtext,
               here,
               janitor,
               metafor,
               pander,
               patchwork,
               tidyverse)

options(DT.options = list(rownames = FALSE,
                          dom = "Blfrtip",
                          scrollX = TRUE,
                          pageLength = 5,
                          columnDefs = list(list(targets = '_all', 
                                                 className = 'dt-center')),
                          buttons = c('copy', 'csv', 'excel', 'pdf')))

source("layout.R")

Custom functions

We also provide some additional helper functions to calculate effect sizes, process data, and visualise our results. The most straightforward way to use these custom functions is to run the code chunk below. Alternatively, paste the code into the console and hit Enter to have R ‘learn’ these custom functions.

If you want to use these custom functions in your own data, you will need to change the variable names according to your own data (check out the R code and you will see what we mean).

# calculate effect sizes ----
## skewness ----
calc.skewness <- function(x, output = "est") {
  n <- length(x)
  
  if (output == "est") { # skewness estimate
    (sqrt(n * (n - 1)) / (n - 2)) *
      (((1 / n) * sum((x - mean(x)) ^ 3)) /
         (((1 / n) * sum((x - mean(x)) ^ 2)) ^ (3/2)))
    
  } else if (output == "var") { # skewness sampling variance
    (6 * n * (n - 1)) /
      ((n - 2) * (n + 1) * (n + 3))
  }
}

## kurtosis ----
calc.kurtosis <- function(x, output = "est") {
  n <- length(x)
  
  if (output == "est") { # kurtosis estimate
    ((((n + 1) * n * (n - 1)) / ((n - 2) * (n - 3))) *
       (sum((x - mean(x)) ^ 4) / (sum((x - mean(x)) ^ 2) ^ 2))) -
      (3 * ((n - 1) ^ 2) / ((n - 2) * (n - 3)))
  } else if (output == "var") { # kurtosis sampling variance
    (24 * n * ((n - 1) ^ 2)) /
      ((n - 3) * (n - 2) * (n + 3) * (n + 5))
  }
}

## Zr ----
r.to.zr <- # Zr estimate
  function(r) { 
    0.5 * log((1 + r) / (1 - r))
  }

zr.variance <- # Zr variance 
  function(n) {
    1 / (n - 3)
  }

## other effect sizes (lnRR and lnVR) ----
calc.effect <- function(data = raw_data, 
                        m) { # calculates other already established effect size statistics 
  escalc(measure = m,
         m1i = data$mean_male,
         m2i = data$mean_female,
         sd1i = data$sd_male,
         sd2i = data$sd_female,
         n1i = data$n_male,
         n2i = data$n_female,
         var.names = c(paste0(m,
                              "_est"),
                       paste0(m,
                              "_var")))
}

# processing functions ----
process.ind_effects <- function(chosen_trait = "fat_mass",
                                measure = "KU_delta") {
  ind_effects <-
    df_meta_analysed %>% 
    filter(trait_name == chosen_trait,
           phenotyping_center %in% c("CCP-IMG",
                                     "HMGU",
                                     "JAX",
                                     "MRC H",
                                     "TCP")) %>% 
    mutate(type = "individual") %>% 
    select(phenotyping_center,
           strain_fig,
           n = n_total,
           est = paste0(measure, "_", "est"),
           var = paste0(measure, "_", "var"),
           lower = paste0(measure, "_", "lower"),
           upper = paste0(measure, "_", "upper"))
  
  model <- rma.mv(data = ind_effects,
                  yi = est,
                  V = var,
                  test = "t",
                  random = list(~ 1|phenotyping_center, 
                                ~ 1|strain_fig))
  
  df_model <- data.frame(trait_name = chosen_trait,
                         est = model$beta[1],
                         var = model$se ^ 2,
                         lower = model$ci.lb,
                         upper = model$ci.ub,
                         phenotyping_center = "Mean",
                         strain_fig = "ES")
  
  
  ind_effects %>% 
    bind_rows(df_model) %>% 
    mutate(est_type = measure,
           centre_and_strain = factor(paste0(phenotyping_center,
                                             "\n",
                                             strain_fig))) %>% 
    mutate(centre_and_strain = factor(centre_and_strain,
                                      levels = c("Mean\nES",
                                                 rev(levels(centre_and_strain)[-6]))))
}

process.cor_effects <- function(chosen_trait_1 = "fat_mass",
                                chosen_trait_2 = "heart_weight") {
  df_effects_cor <-
    df_raw %>% 
    filter(trait_name %in% c(chosen_trait_1,
                             chosen_trait_2), 
           phenotyping_center %in% c("CCP-IMG",
                                     "HMGU",
                                     "JAX",
                                     "MRC H",
                                     "TCP")) %>% 
    pivot_wider(id_cols = c(specimen_id,
                            strain_fig,
                            phenotyping_center,
                            sex),
                names_from = trait_name) %>%
    clean_names() %>% 
    drop_na() %>% 
    group_by(strain_fig,
             phenotyping_center,
             sex) %>% 
    group_modify(~ correlate(.x)) %>% 
    drop_na(all_of(chosen_trait_2)) %>% 
    ungroup() %>%
    left_join(df_raw %>% 
                filter(trait_name %in% c(chosen_trait_1,
                                         chosen_trait_2), 
                       phenotyping_center %in% c("CCP-IMG",
                                                 "HMGU",
                                                 "JAX",
                                                 "MRC H",
                                                 "TCP")) %>% 
                pivot_wider(id_cols = c(specimen_id,
                                        strain_fig,
                                        phenotyping_center,
                                        sex),
                            names_from = trait_name) %>%
                clean_names() %>% 
                drop_na() %>% 
                group_by(strain_fig,
                         phenotyping_center,
                         sex) %>% 
                summarise(n = n())) %>% 
    rename(r_est = chosen_trait_2) %>% 
    mutate(zr_est = r.to.zr(r_est),
           zr_var = zr.variance(n)) %>% 
    select(- c(4:6)) %>% 
    pivot_wider(names_from = sex,
                values_from = c(n,
                                zr_est,
                                zr_var)) %>% 
    mutate(delta_zr_est = zr_est_male - zr_est_female,
           delta_zr_var = zr_var_male + zr_var_female,
           delta_zr_upper = delta_zr_est + 
             qt(0.975, n_male + n_female - 2) * 
             sqrt(delta_zr_var),
           delta_zr_lower = delta_zr_est - 
             qt(0.975, n_male + n_female - 2) *
             sqrt(delta_zr_var))
  
  mlma_zr <-
    rma.mv(data = df_effects_cor,
           yi = delta_zr_est,
           V = delta_zr_var,
           test = "t",
           random = list(~ 1|phenotyping_center, 
                         ~ 1|strain_fig))
  
  df_model <- data.frame(delta_zr_est = mlma_zr$beta[1],
                         delta_zr_lower = mlma_zr$ci.lb,
                         delta_zr_upper = mlma_zr$ci.ub,
                         phenotyping_center = "Mean",
                         strain_fig = "ES")
  
  df_effects_cor %>% 
    bind_rows(df_model) %>% 
    mutate(centre_and_strain = factor(paste0(phenotyping_center,
                                             "\n",
                                             strain_fig))) %>% 
    mutate(centre_and_strain = factor(centre_and_strain,
                                      levels = c("Mean\nES",
                                                 rev(levels(centre_and_strain)[-5]))))
}

# visualisation functions ----
caterpillar.custom <- 
  function(chosen_trait = "fat_mass",
           measure = "KU_delta") {
    plot <-
      process.ind_effects(chosen_trait = chosen_trait,
                          measure = measure) %>% 
      ggplot(aes(y = centre_and_strain,
                 x = est,
                 xmax = upper,
                 xmin = lower,
                 shape = strain_fig,
                 col = phenotyping_center)) +
      geom_pointrange() +
      geom_vline(xintercept = 0,
                 linetype = "dotted") +
      theme_classic() +
      theme(legend.position = "none",
            axis.text.y = element_blank(),
            axis.title.y = element_blank(),
            plot.tag.position = c(0.15, 0.98))
    
    if (measure == "ROM") {
      plot +
        labs(x = "lnRR") +
        scale_x_continuous(limits = c(-0.51, 0.51),
                           breaks = c(-0.5, 0, 0.5)) +
        theme(axis.title.x = ggtext::element_markdown(face = "italic"))
    } else if (measure == "VR") {
      plot +
        labs(x = "lnVR") +
        scale_x_continuous(limits = c(-1, 1),
                           breaks = c(-1, 0, 1)) +
        theme(axis.title.x = ggtext::element_markdown(face = "italic"))
    } else if (measure == "SK_delta") {
      plot +
        labs(x = "&Delta;*sk*") +
        scale_x_continuous(limits = c(-2.1, 2.1),
                           breaks = c(-2, 0, 2)) +
        theme(axis.title.x = ggtext::element_markdown())
    } else if (measure == "KU_delta") {
      plot +
        labs(x = "&Delta;*ku*") +
        scale_x_continuous(limits = c(-15, 15),
                           breaks = c(-15, 0, 15)) +
        theme(axis.title.x = ggtext::element_markdown())
    }
  }

ridgeline.custom <- function(chosen_trait = "fat_mass") {
  df_raw %>% 
    filter(trait_name == chosen_trait,
           phenotyping_center %in% c("CCP-IMG",
                                     "HMGU",
                                     "JAX",
                                     "MRC H",
                                     "TCP")) %>% 
    add_row(phenotyping_center = "Mean",
            strain_fig = "ES") %>% 
    mutate(centre_and_strain = factor(paste0(phenotyping_center,
                                             "\n",
                                             strain_fig))) %>% 
    mutate(centre_and_strain = factor(centre_and_strain,
                                      levels = c("Mean\nES",
                                                 rev(levels(centre_and_strain)[-5]))),
           value_s = scale(value)) %>%
    ggplot(aes(x = value_s,
               y = centre_and_strain,
               fill = sex,
               linetype = sex)) +
    stat_slab(scale = 0.7, 
              alpha = 0.4,
              linewidth = 0.6,
              col = "black") +
    scale_fill_manual(values = c("white",
                                 "black")) +
    scale_linetype_manual(values = c("solid",
                                     "dashed")) +
    labs(x = paste0(str_to_sentence(str_replace_all(chosen_trait,
                                                    "_",
                                                    " ")),
                    "\n(scaled)"),
         y = "Phenotyping centre and mice strain") +
    theme_classic() +
    theme(legend.position = "none",
          axis.title.x = element_text(size = 12, 
                                      margin = margin(t = 0.2,
                                                      unit = "cm")),
          axis.title.y = element_text(size = 12, 
                                      margin = margin(r = 0.2,
                                                      unit = "cm")),
          axis.text.x = element_text(size = 10),
          axis.text.y = element_text(size = 10),
          plot.tag.position = c(0.53, 0.98))
}

cor.caterpillar.custom <- 
  function(chosen_trait_1 = "fat_mass",
           chosen_trait_2 = "heart_weight") {
    
    process.cor_effects(chosen_trait_1 = chosen_trait_1,
                        chosen_trait_2 = chosen_trait_2) %>% 
      ggplot(aes(y = centre_and_strain,
                 x = delta_zr_est,
                 xmax = delta_zr_upper,
                 xmin = delta_zr_lower,
                 shape = strain_fig,
                 col = phenotyping_center)) +
      geom_pointrange() +
      geom_vline(xintercept = 0,
                 linetype = "dotted") +
      labs(y = "Phenotyping centre and mice strain",
           x = "&Delta;*Zr*", 
           shape = "Strain") +
      scale_x_continuous(limits = c(-1, 1),
                         breaks = c(-1, 0, 1)) +
      theme_classic() +
      theme(legend.position = "none",
            axis.title.x = ggtext::element_markdown(size = 12, 
                                                    margin = margin(t = 0.2,
                                                                    unit = "cm")),
            axis.title.y = element_text(size = 12,
                                        margin = margin(r = - 0.1,
                                                        unit = "cm")),
            axis.text.x = element_text(size = 10),
            axis.text.y = element_text(size = 10),
            plot.tag.position = c(0.3, 0.99))
  }

cor.plot.custom <- 
  function(chosen_trait_1 = "fat_mass",
           chosen_trait_2 = "heart_weight",
           chosen_lims = c(-3, 5)) {
    df_cor <-
      df_raw %>% 
      filter(trait_name %in% c(chosen_trait_1,
                               chosen_trait_2), 
             phenotyping_center %in% c("CCP-IMG",
                                       "HMGU",
                                       "JAX",
                                       "MRC H",
                                       "TCP")) %>% 
      pivot_wider(id_cols = c(specimen_id,
                              strain_fig,
                              phenotyping_center,
                              sex),
                  names_from = trait_name) %>%
      clean_names() %>% 
      drop_na() %>% 
      mutate(centre_and_strain = factor(paste0(phenotyping_center,
                                               strain_fig))) %>% 
      mutate(centre_and_strain = factor(centre_and_strain,
                                        levels = rev(levels(centre_and_strain))),
             trait_1_s = scale(get(chosen_trait_1))[,1],
             trait_2_s = scale(get(chosen_trait_2))[,1])
    
    plot_list <- list()
    
    for (i in 1:length(levels(df_cor$centre_and_strain))) {
      level_i <- sort(levels(df_cor$centre_and_strain))[i]
      
      plot <-
        df_cor %>% 
        filter(centre_and_strain == level_i) %>% 
        ggplot(aes(x = trait_1_s,
                   y = trait_2_s,
                   shape = sex,
                   linetype = sex)) +
        geom_point(
          alpha = 0.008,
        ) +
        geom_abline(intercept = 0,
                    slope = 1,
                    linewidth = 0.5,
                    linetype = "dotted") +
        geom_smooth(method = "lm",
                    se = F,
                    col = "black") +
        scale_shape_manual(values = c(3, 4)) +
        scale_linetype_manual(values = c("solid",
                                         "dashed")) +
        scale_x_continuous(limits = chosen_lims) +
        scale_y_continuous(limits = chosen_lims) +
        labs(x = paste0(str_to_sentence(str_replace_all(chosen_trait_1, 
                                                        "_", 
                                                        " ")),
                        "\n(scaled)"),
             y = paste0(str_to_sentence(str_replace_all(chosen_trait_2, 
                                                        "_", 
                                                        " ")),
                        " (scaled)")) +
        theme_classic() +
        theme(legend.position = "none",
              plot.tag.position = c(0.05, 0.91),
              axis.title.x = element_text(size = 12, 
                                          margin = margin(t = 0.2,
                                                          unit = "cm")),
              axis.title.y = element_text(size = 12, 
                                          margin = margin(r = 0.2,
                                                          unit = "cm")),
              axis.text.x = element_text(size = 10),
              axis.text.y = element_text(size = 10))
      
      
      if (i != 6) {
        plot <-
          plot +
          theme(axis.title.x = element_blank(),
                axis.text.x = element_blank(),
                axis.line.x = element_blank(),
                axis.ticks.x = element_blank())
      }
      
      plot_list[[i]] <- plot
    }
    
    return(plot_list)
  }

Equations and custom functions to calculate effect sizes

Skewness

Following Pick et al. (2022).

\[ sk = \frac{\frac{1}{n} \sum_{i = 1}^{n}(x_{i} - \bar{x}) ^ 3}{[\frac{1}{n} \sum_{i = 1}^{n}(x - \bar{x}) ^ 2] ^ \frac{3}{2}} \frac{\sqrt{n (n - 1)}}{n - 2} \] \[ s^2_{sk} = \frac{6n(n - 1)}{(n - 2)(n + 1)(n + 3)} \]

\[ \Delta sk = sk_{1} - sk_{2} \]

\[ s^2_{\Delta sk} = s^2_{sk_1} + s^2_{sk_2} - 2 \rho_{sk} s_{sk_1} s_{sk_2} \]

Kurtosis

\[ ku = \frac{n (n + 1) (n - 1)}{(n - 2)(n - 3)} \frac{\sum_{i = 1}^{n}(x_{i} - \bar{x}) ^ 4} {[\sum_{i = 1}^{n}(x_{i} - \bar{x}) ^ 2]^ 2} - \frac{3(n - 1) ^ 2}{(n - 2)(n - 3)} \] \[ s^2_{ku} = \frac{24 n (n - 1) ^ 2}{(n - 3)(n - 2)(n + 3)(n + 5)} \]

\[ \Delta ku = ku_{1} - ku_{2} \]

\[ s^2_{\Delta ku} = s^2_{ku_1} + s^2_{ku_2} - 2 \rho_{ku} s_{ku_1} s_{ku_2} \]

Zr

\[ Zr = \frac{ln(\frac{1 + r}{1 - r})}{2} \]

\[ s^2_{Zr} = \frac{1}{n - 3} \] \[ \Delta Zr = Zr_{1} - Zr_{2} \]

\[ s^2_{\Delta Zr} = s^2_{Zr_1} + s^2_{Zr_2} -2 \rho_{Zr} s_{Zr_1} s_{Zr_2} \]

Data loading and preparation

We use data from the International Mouse Phenotyping Consortium (IMPC, version 18.0; Dickinson et al., 2016; http://www.mousephenotype.org/).

# raw data ----
df_raw <- 
  read_csv("mice_data_sample.csv") %>% 
  # small adjustments to make plots more readable:
  mutate(phenotyping_center = 
           ifelse(phenotyping_center == "MRC Harwell",
                  "MRC H",
                  phenotyping_center),
         strain_fig = case_when(strain_accession_id == "MGI:2159965" ~ 
                                  "N",
                                strain_accession_id == "MGI:2683688" ~ 
                                  "NCrl",
                                strain_accession_id == "MGI:2164831" ~ 
                                  "NTac",
                                strain_accession_id == "MGI:3056279" ~ 
                                  "NJ",
                                strain_accession_id == "MGI:2160139" ~ 
                                  "NJcl"))

df_meta_analysed <-
  df_raw %>% 
  group_by(sex,
           trait_name,
           phenotyping_center,
           strain_fig) %>% 
  summarize(mean = mean(value,
                        na.rm = T),
            sd = sd(value,
                    na.rm = T),
            n = n(),
            SK_est = calc.skewness(value),
            SK_var = calc.skewness(value, output = "var"),
            KU_est = calc.kurtosis(value),
            KU_var = calc.kurtosis(value, output = "var")) %>% 
  pivot_wider(id_cols = c(trait_name,
                          phenotyping_center,
                          strain_fig),
              names_from = sex,
              values_from = c(mean:KU_var)) %>% 
  mutate(SK_delta_est = SK_est_male - SK_est_female,
         SK_delta_var = SK_var_male + SK_var_female,
         KU_delta_est = KU_est_male - KU_est_female,
         KU_delta_var = KU_var_male + KU_var_female) %>% 
  bind_cols(calc.effect(., m = "ROM")) %>% # lnRR
  bind_cols(calc.effect(., m = "CVR")) %>% # lnCVR
  bind_cols(calc.effect(., m = "VR")) %>% # lnVR
  filter(!is.na(CVR_est)) %>% 
  mutate(n_total = n_female + n_male,
         prop_females = n_female / (n_female + n_male)) %>% 
  select(trait_name,
         phenotyping_center,
         strain_fig,
         n_total,
         prop_females,
         ROM_est,
         ROM_var,
         CVR_est,
         CVR_var,
         VR_est,
         VR_var,
         SK_delta_est,
         SK_delta_var,
         KU_delta_est,
         KU_delta_var) %>% 
  mutate(ROM_upper = ROM_est + qt(0.975, 
                                  n_total - 1) * sqrt(ROM_var),
         ROM_lower = ROM_est - qt(0.975, 
                                  n_total - 1) * sqrt(ROM_var),
         CVR_upper = CVR_est + qt(0.975, 
                                  n_total - 1) * sqrt(CVR_var),
         CVR_lower = CVR_est - qt(0.975, 
                                  n_total - 1) * sqrt(CVR_var),
         VR_upper = VR_est + qt(0.975, 
                                n_total - 1) * sqrt(VR_var),
         VR_lower = VR_est - qt(0.975, 
                                n_total - 1) * sqrt(VR_var),
         SK_delta_upper = SK_delta_est + qt(0.975, 
                                            n_total - 1) * sqrt(SK_delta_var),
         SK_delta_lower = SK_delta_est - qt(0.975, 
                                            n_total - 1) * sqrt(SK_delta_var),
         KU_delta_upper = KU_delta_est + qt(0.975, 
                                            n_total - 1) * sqrt(KU_delta_var),
         KU_delta_lower = KU_delta_est - qt(0.975, 
                                            n_total - 1) * sqrt(KU_delta_var))

Meta-analytical models

We then use the data from multiple phenotyping centres and mice strains to calculate average effect sizes (\(\Delta sk\), \(\Delta ku\), and \(\Delta Zr\)).

Single variable effect sizes

map2_dfr(.x = rep(c("fat_mass",
                    "heart_weight",
                    "glucose",
                    "total_cholesterol"),
                  each = 4),
         .y = rep(c("ROM",
                    "VR",
                    "SK_delta",
                    "KU_delta"), 
                  4),
         .f = process.ind_effects) %>% 
  mutate(est_type = case_when(est_type == "ROM" ~ "lnRR",
                              est_type == "VR" ~ "lnVR",
                              est_type == "SK_delta" ~ "delta_sk",
                              est_type == "KU_delta" ~ "delta_ku")) %>% 
  datatable(.,
            extensions = "Buttons",
            rownames = FALSE)

Correlational effect sizes

map2_dfr(.x = c("fat_mass",
                "glucose"),
         .y = c("heart_weight",
                "total_cholesterol"),
         .f = process.cor_effects) %>% 
  mutate(relationship = rep(c("fat mass and heart weight",
                              "glucose and total cholesterol"),
                            each = 7)) %>% 
  relocate(relationship) %>%  
  datatable(.,
            extensions = "Buttons",
            rownames = FALSE)
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(chosen_trait_2)
## 
##   # Now:
##   data %>% select(all_of(chosen_trait_2))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Visualisations

## figure_2 ----
list_figure_2 <- list()
list_figure_2[[1]] <- ridgeline.custom("fat_mass")
list_figure_2[2:5] <- map2(.x = rep("fat_mass", 4),
                           .y = c("ROM",
                                  "VR",
                                  "SK_delta",
                                  "KU_delta"),
                           .f = caterpillar.custom)
list_figure_2[[6]] <- ridgeline.custom("heart_weight")
list_figure_2[7:10] <- map2(.x = rep("heart_weight", 4),
                            .y = c("ROM",
                                   "VR",
                                   "SK_delta",
                                   "KU_delta"),
                            .f = caterpillar.custom)

(figure_2 <-
    list_figure_2 %>% 
    wrap_plots(ncol = 5) +
    plot_annotation(tag_levels = "A"))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`stat_slabinterval()`).
## Removed 1 row containing missing values or values outside the scale range
## (`stat_slabinterval()`).


## figure_3 ----
list_figure_3 <- list()
list_figure_3[1:6] <- cor.plot.custom(chosen_trait_1 = "fat_mass",
                                      chosen_trait_2 = "heart_weight")
list_figure_3[[7]] <- cor.caterpillar.custom(chosen_trait_1 = "fat_mass",
                                             chosen_trait_2 = "heart_weight")

list_figure_3[8:13] <- cor.plot.custom(chosen_trait_1 = "glucose",
                                       chosen_trait_2 = "total_cholesterol")

list_figure_3[[14]] <- cor.caterpillar.custom(chosen_trait_1 = "glucose",
                                              chosen_trait_2 = "total_cholesterol")
(figure_3 <-
    list_figure_3 %>% 
    wrap_plots() +
    plot_layout(design = layout_2,
                heights = c(rep(1, 6), 0.6),
                widths = c(rep(0.23, 2), 0.02, rep(0.23, 2)),
                axes = "collect",
                guides = "collect") +
    plot_annotation(tag_levels = list(c("A", 
                                        rep("", 
                                            5), 
                                        "B",
                                        "C",
                                        rep("", 
                                            5),
                                        "D"))))
## Warning: Removed 4 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 27 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 27 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 25 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 25 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).


## figure_4 ----
list_figure_4 <- list()
list_figure_4[[1]] <- ridgeline.custom("glucose")
list_figure_4[2:5] <- map2(.x = rep("glucose", 4),
                           .y = c("ROM",
                                  "VR",
                                  "SK_delta",
                                  "KU_delta"),
                           .f = caterpillar.custom)
list_figure_4[[6]] <- ridgeline.custom("total_cholesterol")
list_figure_4[7:10] <- map2(.x = rep("total_cholesterol", 4),
                            .y = c("ROM",
                                   "VR",
                                   "SK_delta",
                                   "KU_delta"),
                            .f = caterpillar.custom)

(figure_4 <-
    list_figure_4 %>% 
    wrap_plots(ncol = 5) +
    plot_annotation(tag_levels = "A"))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`stat_slabinterval()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`stat_slabinterval()`).

Software and package versions

sessionInfo() %>% 
  pander()

R version 4.5.1 (2025-06-13)

Platform: aarch64-apple-darwin20

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: lubridate(v.1.9.4), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.4), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.3.0), ggplot2(v.3.5.2), tidyverse(v.2.0.0), patchwork(v.1.3.1), pander(v.0.6.6), metafor(v.4.8-0), numDeriv(v.2016.8-1.1), metadat(v.1.4-0), Matrix(v.1.7-3), janitor(v.2.2.1), here(v.1.0.1), ggtext(v.0.1.2), ggdist(v.3.3.3), DT(v.0.33), corrr(v.0.4.4) and pacman(v.0.5.1)

loaded via a namespace (and not attached): gtable(v.0.3.6), xfun(v.0.52), bslib(v.0.9.0), htmlwidgets(v.1.6.4), lattice(v.0.22-7), mathjaxr(v.1.8-0), tzdb(v.0.5.0), vctrs(v.0.6.5), tools(v.4.5.1), crosstalk(v.1.2.1), generics(v.0.1.4), parallel(v.4.5.1), pkgconfig(v.2.0.3), RColorBrewer(v.1.1-3), distributional(v.0.5.0), lifecycle(v.1.0.4), compiler(v.4.5.1), farver(v.2.1.2), snakecase(v.0.11.1), litedown(v.0.7), htmltools(v.0.5.8.1), sass(v.0.4.10), yaml(v.2.3.10), pillar(v.1.10.2), crayon(v.1.5.3), jquerylib(v.0.1.4), cachem(v.1.1.0), nlme(v.3.1-168), commonmark(v.1.9.5), tidyselect(v.1.2.1), digest(v.0.6.37), stringi(v.1.8.7), splines(v.4.5.1), labeling(v.0.4.3), rprojroot(v.2.0.4), fastmap(v.1.2.0), grid(v.4.5.1), cli(v.3.6.5), magrittr(v.2.0.3), withr(v.3.0.2), scales(v.1.4.0), bit64(v.4.6.0-1), timechange(v.0.3.0), rmarkdown(v.2.29), bit(v.4.6.0), hms(v.1.1.3), evaluate(v.1.0.4), knitr(v.1.50), mgcv(v.1.9-3), markdown(v.2.0), rlang(v.1.1.6), gridtext(v.0.1.5), Rcpp(v.1.0.14), glue(v.1.8.0), xml2(v.1.3.8), vroom(v.1.6.5), jsonlite(v.2.0.0) and R6(v.2.6.1)

References

Komsta, Lukasz, and Frederick Novomestky. 2022. Moments: Moments, Cumulants, Skewness, Kurtosis and Related Tests. https://doi.org/10.32614/CRAN.package.moments.